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Abstract — We consider a general swarm model of self- 
propelling agents interacting through a pairwise potential in 
the presence of noise and communication time delay. Previous 
work [Phys. Rev. E 77, 035203(R) (2008)] has shown that a 
communication time delay in the swarm induces a pattern 
bifurcation that depends on the size of the coupling amplitude. 
We extend these results by completely unfolding the bifurca- 
tion structure of the mean field approximation. Our analysis 
reveals a direct correspondence between the different dynamical 
behaviors found in different regions of the coupling-time delay 
plane with the different classes of simulated coherent swarm 
patterns. We derive the spatio-temporal scales of the swarm 
structures, and also demonstrate how the complicated interplay 
of coupling strength, time delay, noise intensity, and choice of 
initial conditions can affect the swarm. In particular, our studies 
show that for sufficiently large values of the coupling strength 
and/or the time delay, there is a noise intensity threshold that 
forces a transition of the swarm from a misaligned state into 
an aligned state. We show that this alignment transition exhibits 
hysteresis when the noise intensity is taken to be time dependent. 

Index Terms — Autonomous agents. Delay systems. Pattern 
formation, Nonlinear dynamical systems, Bifurcation 

I. Introduction 

The rich dynamic behavior of interacting multi-agent, or 
particle, systems has been the focus of numerous recent 
studies. These multi-particle systems are capable of self- 
organization, as shown by the various coherent conformations 
with complex structure that they generate, even when the in- 
teractions are short range and in the absence of a leader agent. 
The study of these 'swarming' or 'herding' systems has had 
many interesting biological applications which have resulted in 
a better understanding of the spatio-temporal patterns formed 
by bacterial colonies, fish, birds, locusts, ants, pedestrians, etc. 
IUJ-ItI]. The mathematical study of these swarming systems is 
also helpful in the understanding of oscillator synchronization, 
as in the neural phenomenon of central pattern generators . 
The results of these studies have impacted and have been 
successfully applied in the design of systems of autonomous, 
inter-communicating robotic systems ll9l-ll2ll. as well as mobile 
sensor networks 11311 . 
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It is possible to design swarming models for robotic mo- 
tion planning, consensus and cooperative control, and spatio- 
temporal formation. Pairwise potentials for individual agents 
can be straightforwardly ported onto autonomous vehicles. 
Furthermore, these pairwise interactions can be used in con- 
junction with simple scalable algorithms to achieve multi- 
vehicle cooperative motion ll4ll . Specific goals include: ob- 
stacle avoidance lllll . boundary tracking Iil5il , environmental 
sensing 11131 11611 and decentralized target tracking lIlTIl . 

An important problem is that of environmental consensus 
estimation. Here, the individuals of the swarm communicate 
with each other through a network to achieve asymptoti- 
cally synchronous information about their environment 1113 1. 
Recently, consensus was extended to include time delayed 
communication among agents lIlSll . 

Task allocation is another problem of interest involving 
robotic swarms. The objective is to reallocate swarm robots 
to perform a set of tasks in parallel and independently of one 
another in an optimal way. In order to make task reallocation 
more realistic it is possible to consider a time delay that arises 
from the amount of time required to switch between tasks |119 1 . 

Regardless of the design objective of a robotic swarm 
system, a comprehensive theoretical analysis of the model 
must be performed in order to achieve successful algorithm 
design. 

Many different mathematical approaches have been utilized 
to study aggregating agent systems. Some of these studies 
have treated the problem at a single-individual level, using 
ordinary differential equations (ODEs) or delay differential 
equations (DDEs) to describe their trajectories ifiol 120-221. 
An alternative method has been proposed by other researchers 
and consists of using continuum models that consider averaged 
velocity and agent density fields that satisfy partial differential 
equations (PDEs) 12. S S 0] ■ In addition, authors also have 
studied the effects of noise on the swarm's behavior and have 
shown the existence of noise-induced transitions lEsl 251. 



The study of these systems has been enriched by tools from 
statistical physics since both first and second order phase 
transitions have been found in the formation of coherent states 



An additional effect that has recently been considered is 
that of communication time delays between robots. Time delay 
models are common in many areas of mathematical biology 
including population dynamics, neural networks, blood cell 
maturation, virus dynamics and genetic networks 11291-1371. 
In the context of swarming agents, it has been shown that 
the introduction of a communication time delay may induce 
transitions between different coherent states in a manner which 
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depends on the coupling strength between agents and the noise 
intensity 12511 . Thus far, most of the work has concentrated 
on the case of uniform time delays among agents ll26ll . 
However, the practical engineering of multi-agent systems 
requires researchers to consider the case in which time delays 
may vary due to data processing times, problems in inter-agent 
communication, etc. The case of differing (and even time-vary 
ing) time delays between agents may be treated similarly to 
the case of a single delay by using a data buffer iIztIi . 

In this work, we carry out a detailed study of the bifurcation 
structure of the mean field approximation used in ll25ll and 
investigate how the bifurcations in the system are modified 
in the presence of noise. Section HI] contains the swarm 
model, while Sec. |lll] contains the derivation of the mean 
field approximation. The bifurcation analysis of the mean 
field equation can be found in Sec. IIVI and Sec. IV] provides 
a comparison of the mean field analysis with the nonlinear 
governing equations. In Sec. IVII we describe the effects of 
noise on the swarm, and the conclusions are contained in 

secrvm 

II. Swarm Model 

We consider a two-dimensional (2D) swarm that consists 
of N identical self-propelling individuals of unit mass that 
are mutually attracted to one another in a symmetric fash- 
ion. Hence, the coupling of the agents occurs via a fully 
connected graph. In addition, we consider the case in which 
the individuals that comprise the swarm are communicating 
with each other in a stochastic environment. Because of the 
finite communication times between individuals, there is a time 
delay between interactions. Assuming that the communication 
time between agents is constant and equal to r > 0, the swarm 
dynamics is described by the following governing equations: 



(la) 



N 



= (1 - |v,n V, - - ^(r,(t) - r,(t - r)) + q^, 

(lb) 



for i = 1,2 . . . , N. The terms and respectively represent 
the 2D position and velocity of the i-th agent at time t. 
The strength of the attraction is measured by the coupling 
constant a > 0. The self-propulsion and frictional drag forces 
on each agent is given by the term (l — |vip) v^. Therefore, 
in the absence of coupling, agents tend to move on a straight 
line with unit speed |vi| = 1 as time goes to infinity. The 
term ri^{t) = (f?,-^', f/P') is a 2D vector of stochastic white 
noise with intensity equal to D such that {i]\^\t)) = and 
{v-'\t)rij''\t')) = 2D6{t-t')Sij5ik for i,j = 1, 2, . . . and 
£,k = 1,2. It is the main objective of this work to identify 
the possible swarm behaviors for different values of a and t. 

The coupling between individuals arises from a time de- 
layed, spring-like potential. Hence, our equations of motion 
may be considered to be the first term in a Taylor expansion 
of other more general time delayed potential functions about 
an equilibrium point. 



III. Mean Field Approximation 

We can investigate the stability of the swarm system by 
deriving a mean field approximation of the system. The deriva- 
tion involves the consideration of agent coordinates relative to 
the center of mass and the elimination of the noise terms. The 
center of mass of the swarming system is given by 



1 ^ 



1=1 



The position of each individual can be decomposed into 
r, = R + (5r,, i^l,2...,N, 



(2) 



(3) 



where Svi is the vector from the center of mass to particle i 
and 



N 



(4) 



We substitute the ansatz given by Eq. ([3]) into the second order 
system that is equivalent to Eqs. (fTali-(fTb]i with D = Q. After 
simplification, one obtains 

R + 5r, = (l - |R.|^ - 2R • Sh - \Sri\^^ (R + Svi) 
a{N -I) ( ^^^^ _ ^^^^^^ 



N 



N 



Sr,{t - t), 



(5) 



where we used the fact that Eq. can be written as 



N 



6r,it-T) = -J2Srjit-T). (6) 
Summing Eq. (|5]l over i and using Eq. (|4|i, we find 



R=(^l-|Rp-l^|5f,p^R 



2R-6r, + \6r,\^)5r. 

i=l 

N ~1 



(7) 



By inserting Eq. O into Eq. (|5]l it is possible to find the 
following equation for (Jr^: 



i^f]|<5r,f -2R.5r,-|<5r,pjR 



1 - |R|^ - 2R-(5r, - |(5r,|^j Sr. 

, N 

-^(2R.<5r, + |<5r,p) ^r,- 



N 



N 



N 



(8) 



fori = l,2...,iV. 



Taken together, Eqs. (ITJ and (|8]l are equivalent to Eqs. JTaj - 
( fTbl i and they merely involve a reconstruction of the original 
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system that is written in terms of particle coordinates into 
this new system that is written in terms of the center of mass R 
and coordinates relative to the center of mass 5ri. One can see 
that this mapping has transformed the original 2N differential 
equations into 2N + 2 equations. Due to the relation given 
by Eq. dU, only 2N of the transformed set of equations are 
independent. Therefore, there is no inconsistency between the 
original and transformed equations. 

By neglecting the fluctuation terms 5ri from Eq. (|7J and 
taking — > oo, we obtain the following heuristic mean field 
approximation for the center of mass: 



(9) 



R = 1 - |Rr R - a (R(0 - R(i - t)) , 



where we made the approximation a^^^ ~ a since we are 
considering the large system size limit N ^ oo. We will 
address the validity of neglecting the fluctuation terms in 
Section [V] 



IV. Bifurcations in the Mean Field Equation 

Having derived a mean field equation, we continue by 
analyzing the bifurcation structure. This bifurcation analysis 
will allow us to better understand the behavior of the system 
in different regions of parameter space. Letting R = {X, Y) 
and R — {U, V), Eq. ^ may be written in component form 
as 

X = U, (10a) 
if ^ (l-U^ -V^)U - a{X - X{t-T)), (10b) 

Y = V, (10c) 

V = -V'^)V -a{Y -Y{t-T)). (lOd) 



Regardless of the value of a and r, Eqs. ( llOab -l fTOdt have 
translational invariant stationary solutions given by 



X = Xo, U = Q, Y = Yo, V = 0, 



(11) 



where Xq and Yq are two free parameters. In addition, Eqs. 
( llOab -l fTOdl i also have a three parameter family of uniformly 
translating solutions given by 



X^Uot + Xo, U^Uo, Y^Vot + Yo, V^Vo, 



which requires 



Uq + Fq^ = 1 - ar 



(12) 



(13) 



and is real-valued only when ar < 1. In the two-parameter 
space (a, r), the hyperbola ar = 1 is in fact a pitchfork 
bifurcation curve on which the uniformly translating states 
are bom from the stationary state {Xq, 0, Iq, 0). The pitchfork 
bifurcation curve can be seen in Fig. |l(a)| The other branch 
of the pitchfork bifurcation is an unphysical solution with 
negative speed. 

Linearizing Eqs. jlOab - dlOdb about the stationary state, we 
obtain the characteristic equation 



(a(l-e-^")-A + A2) =0 



(14) 



It is sufficient to study the zeros of the function 

V{X)=a(l-e-^^)-X + X^^O, (15) 

since the eigenvalues [see Eq. ( fT4l il of the system given by 
Eqs. ( llOab -l fTOdl i are obtained by duplicating those of Eq. ( fTSl ). 

We identify the Hopf bifurcations in the two parameter 
space (a, t) by letting the eigenvalue be purely imaginary. 
Our choice of A = ioj is substituted into Eq. ( fTSl l. and one 
obtains 



(16) 



By taking the modulus of Eq. ( fTSl ). one finds that a at the 
Hopf point is given by 



(an 



.2\2 



If we consider the case when uj ^ Q, then 



1 



an = 



(17) 



(18) 



We substitute Eq. ( fTSl l into Eq. ( fT6l ) and take the complex 
conjugate. This allows us to obtain the following equation for 
T at the Hopf point that does not involve a: 



1 



2uj 



j-e'"". (19) 

We isolate r by equating the arguments of both sides, being 
careful to use the branch of tan 6* in (0, tt) since the left hand 
side of the equation above is on the upper complex plane for 
w > 0. We then obtain a family of Hopf bifurcation curves 
parameterized by uj: 



1 



th„ (i^) = — I arctan 

UJ 



(20a) 



2w 



2mT 



0,1,, 



(20b) 



The first few members of the family of Hopf bifurcation curves 
are shown in Fig. |l(a) It also is possible to eliminate the 
parameter ut in Eqs. ( I20al )-( l20bb . Doing so, one obtains 



TH„ (a) = 
1 

V2a- 1 



arctan 



V2a - 1 
1 - a 



2mT 



n = 0,1,. 



(21) 



In spite of their appearance, the Hopf curves in Eqs. 
(I20ab -( l20b] i and (ISTT i are in fact continuous at w = 1 and 
a = 1, respectively [with the correct branch of tan6' in 
(0, tt)]. Inspection of Eq. ( I20ab . shows that the Hopf frequency 
depends only on the value of a for all members in the family. 
The frequency equals one when a = 1, and the frequency tends 
to infinity as a increases. Interestingly, only the first Hopf 
curve is defined at a = 1/2 and has the value TH„\a=i/2 = 2. 
The point (a = 1/2, r = 2) which lies both on the first Hopf 
curve and on the pitchfork curve is a Bogdanov-Takens (BT) 
point (the eigenvalues are zero), where w = 0. None of the 
other Hopf branches meet the pitchfork bifurcation curve since 
r — > oo as a — > 1/2. 
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Fig. 1 . (a) Hopf (blue) and pitchfork (red) bifurcation curves in (a.r) space, 
(b) A zoom-in of Fig. |l(a)| Included is the saddle to node transition curve 
(dashed black) and a number in each region (with boundaries given by the 
solid curves) that indicates the number of eigenvalues with a real part greater 
than zero. 




-0.05 0.05 0.1 

9?(A,) 





-0.05 

9?(X,) 



-0.2 



-0.04 -0.02 0.02 

9?(X) 

Fig. 2. Real and Imaginary parts of the dominating eigenvalues as one moves 
around the Bogdanov-Takens point (a = 1/2, r = 2) in (a, r) parameter 
space. The eigenvalues shown are associated with the locations (a) a = 0.60, 
T = 2.0, (b) a = 0.48, r = 2.09, (c) a = 0.40, r = 2. 01, (d) a = 0.53, 
T = 1.90, and (e) a = 0.55, r = 1.91. Refer to Fig. |l(b)| to see where each 
of the (a, t) points lies in relation to the bifurcation curves. 



The pitchfork and Hopf bifurcation curves in the (a, t) 
parameter space were computed using a numerical continu- 
ation method Issll . These resuhs (not shown) are in perfect 
agreement with our analytical calculations. These numerical 
continuation studies also allow for the determination of the 
number of eigenvalues with real part greater than zero in 
different regions of the (a, t) parameter space. The results are 
shown in Fig. |l(b)| In addition, our numerical continuation 
analysis revealed node to focus transitions of the steady state. 
These transitions occur at points where there are two real and 
equal eigenvalues, i.e. where I?(A) — and I?'(A) = 0, 
for real- valued A. If I?' (A) = then one can show that 
e""^^ = ^^^^ . Insertion of this relation into ViX) = leads 
to 

^0, (22) 



A" 



1 



A + a 



which has solutions A = 



For the 



roots to be repeated, we set the discriminant equal to zero and 
this gives the following curve where the node-focus transitions 
occur: 



1 



(23) 



Moreover, by inspecting the solutions to Eq. (|22] | one finds 
that the repeated eigenvalues have positive real parts if t > 2 



and negative real parts if r < 2. In Figure |l(b)| we show 
the pitchfork and Hopf bifurcation curves overlaid with the 
node-focus transition curve given by Eq. (l23T l. 

As seen in Fig. |l(b)[ the pitchfork and first Hopf bifurcation 
curves, together with the node-focus transition curve, split the 
area around the BT point into five different regions. The be- 
havior of the dominating eigenvalues (excluding the one at the 
origin) in each of these five regions is shown in Figs. |2(a)|2(e)] 
Starting at a point directly to the right of the BT point in (a, t) 
space, there is a pair of eigenvalues with positive real parts 
and non-zero imaginary parts [Fig. |2(a)| . Moving counter- 
clockwise, the eigenvalue pair collapse on the positive real axis 
upon crossing the upper branch of the node-focus transition 
curve [Fig. |2(b)] . Continuing in the same direction, we observe 
two different instances of eigenvalues crossing the origin: (i) 
first the smaller of the two purely real and positive eigenvalues 
does so as the upper part of the pitchfork bifurcation curve is 
crossed [Fig. |2(c)| and (ii) then the remaining purely real and 
positive eigenvalue crosses the origin as the lower part of the 
pitchfork bifurcation curve is crossed [Fig. |2(d)| . Finally, as 
the node-focus transition curve is crossed, the two purely real 
and negative eigenvalues coincide on the negative real axis 
and acquire non-zero imaginary parts [Fig. |2(e)| . Continuing 
upwards in parameter space, the complex pair of eigenvalues 
crosses the imaginary axis as the Hopf bifurcation curve is 
crossed, giving birth to a stable limit cycle. 
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a 

Fig. 3. Regions in (a, t) space with different dynamical behavior. 



V. Comparison of the Mean Field Analysis and the 
Full Swarm Equations 

Our analysis of the deterministic mean field equations iden- 
tified the different dynamical behaviors that the approximation 
given by Eq. ^ exhibits in different regions of the (a, r) 
plane. However, the analysis does not provide any information 
about how the swarm agents are distributed about the center 
of mass. We neglect the stochastic terms in Eqs. (fTali-(frb]i and 
use extensive numerical simulations to identify some of the 
coherent structures that the swarm adopts asymptotically in 
time; 

(i) A translational state, in which all swarm particles have 
identical positions and velocities and move uniformly in 
a straight line. The direction of motion depends on the 
initial conditions. This behavior is only possible in region 
A of Fig. |3] Moreover, the asymptotic convergence to 
this state requires that all particles be located in close 
proximity and with aligned velocities at the initial time. 
Hence, the basin of attraction is extremely small which 
causes this state to be very sensitive to perturbations. This 
is discussed in more detail below. 

(ii) A ring state, in which the center of mass is stationary. The 
swarm agents distribute themselves along the ring with 
roughly half of the agents moving clockwise and half of 
the agents moving counter-clockwise. The final stationary 
position of the center of mass and the particular behavior 
of each individual in the swarm is dependent on the initial 
conditions. This behavior is possible in regions A, B and 
C of Fig. S 

(iii) A rotational state, in which all swarm agents collapse to 
the center of mass and the latter rotates on a circular orbit. 
The direction of rotation depends on the initial conditions. 
This behavior is only possible in region C of Fig. |3] 

(iv) A degenerate rotational state, in which all swarm particles 
collapse to the center of mass and the latter oscillates 
back and forth on a line. This behavior is only possible 
in region C of Fig. |3] In addition, it requires that the 
initial motion of all swarm particles be constrained to a 
line and so is sensitive with respect to perturbations and 
noise. 



The above list is not extensive and our simulations have 
revealed other time-asymptotic patterns. However, all of these 
other patterns (and including the translational state and the 
degenerate rotational state) require extreme symmetry in the 
initial conditions and are very sensitive with respect to pertur- 
bations and noise. Our numerical simulations suggest that only 
the ring and the rotational state have a significant robustness 
with respect to perturbations and noise. 

The full system of equations predict a bistable behavior 
since the translating and ring states are both possible in region 
A and C [Fig. [3], depending on the initial conditions. The 
linear stability analysis of Section |IV] shows that the mean 
field approximation fails to capture this bistable behavior. 

The mean-field bifurcation results obtained here are of prac- 
tical value since they provide us with guidelines for selecting 
values for a and r that will result in a particular coherent 
pattern asymptotically in time. In the case of bistability, 
our numerical simulations strongly suggest that the initial 
alignment of the agents' velocities is critical in determining the 
coherent state adopted. Specifically, to obtain the translating, 
rotating and degenerate rotating states asymptotically in time 
(structures in which the individuals' velocities are perfectly 
aligned), one requires a high alignment of the initial particles' 
velocities; otherwise, the swarm will adopt the ring state. 
However, how high an alignment is needed depends on the 
specific choice of (a, r). Our results indicate that it is easier to 
obtain aligned states for larger values of the coupling constant 
a. Unfortunately, it is not feasible to obtain analytic basin 
boundaries in this infinite dimensional system. In principle, 
one may approximate such boundaries by performing pro- 
hibitively extensive numerical simulations where the space 
of history functions is restricted in some way. Therefore, the 
computation of basins of attraction is outside the scope of this 
work and is left for future research. 

For the non-degenerate and degenerate rotating states as 
well as for the translating state, the approximation we made 
when neglecting the fluctuation terms in Eq. ^ is entirely 
valid since in the noiseless case all agents collapse to the center 
of mass. In the case of the ring structure, these fluctuation 
terms are not necessarily small. However, in Eq. (|7) all 
fluctuation terms with the exception of the one containing the 
factor ;^ X^ili I'^^iP approximately cancel out in the long 
time limit, due to the symmetry in the distribution of the 
agents. The fluctuation term that remains becomes equal to 
one in the long time limit. This has the effect of eliminating 
the self-propulsion of the center of mass and what remains is 
solely cubic dissipation. 

The following sub-sections contain detailed discussion re- 
garding the spatio-temporal scales of each coherent structure. 

A. The Ring State 

The analysis of Appendix |A] shows that the radius and 
angular frequency of the swarm particles on the ring state 
is given by 

Pj = -)=, - ±VH, (24) 

V a 

so that particles move at unit speed, pjOj = ±1. 
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Fig. 4. Comparison of numerical simulations (red circular markers) with the 
analytical expressions (continuous blue curve) given by Eq. )24t for (a) the 
radius and (b) the frequency of the ring state. (c) For each value of a, the time 
delay was chosen as t = a — 1/4 (black circular markers). 




Fig. 5. In (a,r) space, we plot: Hopf (blue) and pitchfork (red) bifurcation 
curves, and the curve ar^ = 2 where the first limit cycle ceases to exist by 
having its radius diverging to infinity (green). 



We have numerically computed the radius and angular 
frequency for different values of a and r within the region in 
which the mean field approximation gives a stable stationary 
center of mass (Fig. |4]i. Figures |4(a)|4(b)| shows that there is 
excellent agreement between the numerical simulations and 
the analytical result given by Eq. (l24l i. It is worth noting that 
the condition given by Eq. ( |29] l and used to derive Eq. (|24] | is 
satisfied in the long time limit in our simulations. 

B. The Rotating State 

We show in AppendixlBlthat the circular orbit of the rotating 
state has radius po and frequency uj that satisfy the following 
relations: 



=a • (1 — cos wr) 
1 



(25a) 

Po =^\/l-a^^^- (25b) 

Eqs. ( |25a| )-( [25b] i can have as many solutions as desired by 
choosing a and r large enough. However, a careful analysis 



Fig. 6. Comparison of numerical simulations (red circular markers) with 
the analytical expressions (continuous blue curve) given by Eqs. j25at - 
)25b) for (a) the radius, (b) the period, and (c) the speed of the collapsed 
circular orbit, (d) For each value of a, the time delay was chosen as 
r = ^ arctan ( ^ ) (black circular markers) to assure asymptotic 

time convergence to tlie collapsed circular orbit state. 



reveals that the solutions to Eqs. (|25at -( [25b] i are generated 
exactly along the Hopf curves of our previous mean field 
analysis and represent the same limit cycles of that analysis 
[Fig. |l(a)| . The expressions in Eqs. (I25ab -( l25bb thus determine 
the spatio-temporal scales of these circular orbits beyond the 
Hopf curves where they are bom. Our analysis also shows that 
the circular limit cycle that is created on the first member of 
the Hopf bifurcation curves persists to the left of the pitchfork 
bifurcation curve and then ceases to exist as its radius diverges 
to infinity on the curve ar^ = 2 (Fig.|5]l. Moreover, numerical 
simulations of the mean field equations reveal that both the 
translating state and the rotating state are linearly stable for 
(a, r) pairs inside the wedge between the curve ar^ — 2 and 
the pitchfork bifurcation curve ar = 1 above the BT point. 

Figures |6(a)|6(d)1 show the excellent agreement between 
numerical simulations and the analytical results given by Eqs. 
(I25al i-( l25b] i. for different values of a and t. 

Interestingly, in Fig. |6(c)| we note that in the asymptotic 
time limit the collapsed agents move at a speed greater than 
one, the speed at which agents would tend to move in the 
absence of coupling. This is explained by noting that the ratio 
of the time delay to the period of oscillations is such that the 
delayed position of the collapsed agents R(i — r) is ahead 
of the present position R(<). The attraction that an individual 
particle feels to the delayed position of the rest of the swarm 
forces the whole system go faster 

C. The Degenerate Rotating State 

A degenerate version of the rotating state is possible when 
the initial motion of the swarm is restricted to a line, since 
in this case it follows from Eqs. (fTali-(fTb]i that the swarm will 
remain on such a line for all times. As we show in Appendix IC] 
we may assume that the motion of the collapsed swarm occurs 
on the X = Y line of the center of mass coordinates and 
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Fig. 7. Comparison of numerical simulation (red circular' markers) with the 
analytical expressions (continuous blue curve) given by Eqs. <42al-<42bl for 
(a) the amplitude, (b) period, and (c) the maximum speed of the collapsed 
straight line orbit, (d) At each value of a, the time delay was chosen as 
(black circular markers) to ensure asymptotic 



/2a-l " V 1- 

trme convergence to the collapsed-straight line orbit state. 



then use a finite Fourier mode approximation of the ensuing 
dynamics. An approximation in terms of just three modes gives 

X{t) = Y{t) = 2ci cosujt + 2|c3| cos(3wt + ^3), (26) 



where uj, c\, C3 and 03 are obtained by solving Eqs. ( 142 ab - 
(I42bl l numerically. 

Figures [7(a)|7(d)| show a comparison between our analytical 
results given by Eqs. ( I42ab - (l42bb and results obtained using 
numerical simulation for the amplitude, period and maximum 
speed of oscillation for different values of a and t. There 
is excellent agreement in both amplitude and period between 
our analysis and the numerical simulations [Figs. |7(a)|7(b)) . 
The agreement for the speed of motion is very good as well, 
but the theoretical estimate is shifted slightly with respect to 
the results from simulations [Fig. |7(c)) . As in the collapsed 
circular orbit, we note that the collapsed set of particles have a 
maximum speed which exceeds one, the speed that individual, 
uncoupled particles acquire in the long-time limit. As before, 
this effect arises from the attraction that the current particle 
position R(i) feels towards the delayed position R(t — r) 
when the latter lies in the direction of motion of the collapsed 
particles. 

VI. The Effects of Noise on the Swarm 

In the absence of noise, the initial alignment of the swarm 
particles is critical in determining the asymptotic behavior of 
the swarm (Sec. |V]i. When noise is introduced, the interplay 
of coupling strength, time delay and noise intensity gives 
rise to very interesting behavior due to fluctuations in the 
particles' alignment. Specifically, our studies show that if the 
coupling strength a and/or the time delay t are below a certain 
limit, then the presence of noise promotes swarm transitions 
from aligned into misaligned coherent states. More surprising, 
however, is that if the coupUng strength a and/or the time delay 



T are big enough, then there is a noise intensity threshold 
that forces a transition in the swarm from misaligned into 
aligned states. In addition, we show that for these high values 
of a and/or r, the system presents an interesting hysteresis 
phenomenon when the noise intensity is time dependent. 

For the purpose of these studies, we define the alignment 
of particle j with the rest of the swarm as the cosine of the 
angle between the velocity of particle j and the velocity of 
the swarm as a whole; 



cos 



f, R 



IRI 



(27) 



Therefore the alignment of individual particles ranges from -1 
to 1. A good measure of the overall alignment of the swarm 
is furnished by the ensemble average of these cosines given 
as 

JV 



Mean swarm aUgnment = — cos 6* , = — 



1 V- r , R 



N 



(28) 

We first carry out a numerical simulation with coupling 
constant a = 0.5 and noise standard deviation a = 0.05 (noise 
intensity D = 0.00125). Al t = 50, a time delay of r = 0.5 
is turned on. These parameters correspond to region A of Fig. 
|3] Initially, we place all particles at the origin and align their 
velocities by choosing Xj = 1 and yj = 1 for all particles. We 
describe the behavior of the swarm by following the ensemble 
averages of the particle distances to the center of mass [Fig. 
|8(a)| and of the particle alignment [Fig. |8(b)| as functions 
of time. Before the time delay is turned on at i = 50, the 
swarm is in a translating state with particles slightly spread 
out from the center of mass in a 'pancake' shape, as described 
in I23II . with an ensemble alignment close to one. Once the 
delay is turned on, the translating state is broken up and the 
swarm converges to the ring state in which the mean particle 
alignment is near zero. The radius of the ring obtained in this 
numerical simulation matches the theoretical result [Eq. d24l l1 
that predicts a radius of = -\/2 « 1.41. A completely 
analogous situation ensues for parameters in region B of Fig. 
[3](results not shown). In addition, in both cases the swarm will 
immediately converge to the ring state if the swarm velocities 
are not sufficiently aligned at time zero. We thus conclude 
that for these choices of (a, r) pairs, the noise misaligns the 
particles' velocities and forces a transition into the ring state. 

In contrast to the cases discussed above, for parameters 
in region C of Fig. [3] a sufficiently large noise intensity 
promotes transitions from misaligned to aligned states. We 
show this by comparing the results of a series of simulations 
for different values of the noise standard deviation cr. The 
simulations are divided into two cases that differ only on the 
initial conditions for the swarm particles. In all simulations, 
the coupling constant a = 2 and a time delay of r = 2 is 
turned on at i = 50. In the first case, all particles start from 
the origin with identical velocities Xj — 1 and yj = 1. In 
the second case, all swarm particles are initially distributed 
uniformly on the unit square and are at rest. 

In these simulations, the final state of the swarm may be vi- 
sualized by plotting the mean swarm alignment after transients 
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Fig. 8. Time evolution of the ensemble average of (a) the particle distance 
to the center of mass, and (b) the mean particle alignment showing how the 
particle alignment breaks up due to the effects of noise. For long times the 
swarm converges to a ring state. The parameter values of a = 0.5 and r = 0.5 
are associated with region A of Fig. |3] The time delay is turned on at t = 50 
and the noise standard deviation is cr = 0.05 {D = 0.00125). 



States of the swarm demonstrate that the system exhibits a 
hysteresis phenomena. With the swarm system starting on the 
ring state with noise standard deviation of a = 0.24, one can 
force a transition into the rotating state by increasing the noise 
to <T = 0.26. However, even if the noise is lowered down to 
(J = 0.02, the swarm remains in the rotating state with a 
high velocity aUgnment [Figs. |1 l(a)|[TT(b)) . Nevertheless, it is 
possible to return the swarm to the ring state if, once in the 
rotating state, the noise is raised to very high amounts (a = 1) 
for a sufficient amount of time and then dropped suddenly to 
a very low value (cr = 0.05). The high noise levels serve to 
completely misalign the particles' velocities and allow them to 
converge to the ring once the noise levels are below a < 0.26. 
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Fig. 9. Asymptotic value of the mean particle alignment for (a) particles 
starting with perfectly aligned velocities at time zero and (b) for particles 
distiibuted uniformly over the unit square and starting from rest for different 
values of the noise standard deviation a. The parameter values of a = 2 and 
T = 2 (turned on at t = 50) are associated with a location in region C of 
Fig.H] 



Fig. 10. Time evolution of the ensemble average of (a) the particle distance 
to the center of mass, and (b) the mean particle alignment showing how 
the swarm transitions from a ling state into a compact, rotational state with 
aUgnment close to one. The parameter values of a = 2 and t = 2 (turned 
on at t = 50) and cr = 0.4 (D = 0.08) are associated with region C of Fig. 
[5] Particles are initially distributed uniformly over the unit square and start 
from rest. 



have decayed (t = 300) as a function of noise intensity for 
the first case [Fig. |9(a)) and the second case [Fig. [9(b)). In the 
first case of simulations, the high initial alignment of particles' 
velocities forces the swarm to converge to a compact rotating 
state independent of noise intensity. However, the rotational 
state is destroyed if the noise standard deviation is bigger 
than a « 0.8 [Fig. |9(a)| . The situation is more interesting 
and complex for the second set of simulations. For low noise 
intensities (a < 0.26) the low initial alignment of the particles 
leads the swarm to converge to a ring state with near zero mean 
alignment [Fig. |9(b)| . A noise standard deviation just beyond 
the threshold of cr w 0.26 displays an interesting effect. As the 
cr ss 0.26 threshold is crossed, the swarm transitions from the 
ring state into the rotating state with high mean alignment. 
An examination of the full simulation data reveals that the 
transition occurs as an increasing group of particles gradually 
becomes aligned and eventually absorbs all the remaining 
particles. A sufficient amount of noise is necessary for this 
transition, since it allows each particles' velocity vector to 
probe many directions until finally enough of them become 
trapped in a 'potential well' of alignment with other particles. 
As with the first case of simulations a noise standard deviation 
bigger than cr « 0.8 breaks up the rotating state. Figure [TO] 
clearly shows the transition from the ring to the compact, 
rotational state through the time evolution of the ensemble 
averages of the particle distances to the center of mass and of 
the mean particle alignment. 

Further studies on the switching behavior between coherent 



VII. Conclusions 

In this work we analyzed the dynamics of a self-propelling 
swarm where individuals interact with a communication time 
delay in the presence of noise. Using a mean field approx- 
imation in the deterministic case, we analytically obtained 
the complete bifurcation picture in the parameter space of 
coupling strength and communication time delay. This analysis 
shows how different combinations of coupling strength and 
time delay induce the swarm to adopt different coherent 
structures asymptotically in time. Our bifurcation studies 
demonstrated the existence of a Bogdanov-Takens point, where 
the stationary center of mass solution has a double zero 
eigenvalue, which is critical in organizing the dynamics of 
the swarm. 

The stable patterns that are possible for this system have 
several applications for autonomous vehicles. More detailed 
applications for each pattern are as follows: (1) the trans- 
lational state may be used for target tracking and group 
transport 111 U 11711 . (2) The ring state should prove useful 
in terrain coverage and regional surveillance 11391 koll . (3) 
The rotating state may be exploited in obstacle avoidance, 
boundary tracking and surveillance lUI, 15, 40 j. In addition, 
we believe all three patterns are applicable to the problem of 
environmental sensing ll3l 11611 . 

In numerical experiments with noise, we showed that the 
interplay of coupling strength, time delay and noise intensity 
may give rise to interesting switching behavior from one 
coherent structure to another We found that if the coupling 
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Fig. 11. Time evolution of (a) mean particle alignment for example 1 , 
(b) noise standai'd deviation for example 1, (c) mean particle alignment for 
example 2, and (d) noise standard deviation for example 2. The results show 
how a time-dependent noise intensity may be used to force swann transitions. 
The parameter values of a = 2 and r = 2 (turned on at t = 10) are associated 
with region C of Fig. |3] Particles are initially distributed uniformly over the 
unit square and start from rest. 



Strength a and/or the time delay r are below a certain limit, 
then the presence of noise induces transitions from states in 
which the alignment of the particles' velocities is high into 
states with low alignment. More surprising, however, is that if 
the coupling strength a and/or the time delay r are big enough, 
then there is a noise intensity threshold that forces a transition 
in the swarm from misaligned into aligned states. In addition, 
by using a time-dependent noise intensity at these high values 
of a and/or r, we show that the system exhibits hysteresis 
since the swarm's transitions are not easily reversible. We 
note that analytical results on the effects of noise on delay- 
coupled swarms are not easy to obtain. Two examples relevant 
to our work are given in ll23[ 12411 . where the authors investigate 
models similar to the one presented here but without time 
delay. 

Realistic application of the model treated here to the motion 
of multi-robot systems requires local repulsion among individ- 
uals to be taken into account. We have simulated the swarm 
model with the addition of a repulsive inter-agent potential of 



exponential form [/^ 



demonstrate (results not shown) that the coherent patterns 
we discussed in this article persist when the characteristic 
repulsion length and repulsion strength between robots 
are small compared to global attraction parameters. Stronger 
repulsion can destabilize the coherent structures. 

Recently, systems with non-uniform time delays have re- 
ceived much attention. For example, the important question 
of synchronization in networks communicating at randomly- 
distributed time delays has been recently investigated ll4lll42ll . 
In practical applications, the case of differing (and even time- 
varying) time delays between agents may be treated similarly 
to the case of a single delay by using a data buffer 12711 . The 
idea is to identify an upper bound to the time delay (Tmax) 
between all agent pairs and then design the agents so that the 



These simulations 



actuation occurs when the data buffer of size T,nax is full. 

As part of our ongoing work, we are extending our in- 
vestigations for the cases in which: (i) the communication 
time delays vary between different pairs of agents; and (ii) 
the communication graph is non-globally coupled. In realistic 
settings, both of these cases may occur due to the effects of 
the spatial distribution of agents such as signal travel times 
and imperfect transmission arising, for example, from complex 
terrain topography or component malfunction. In the case of 
communication delays that differ among different pairs of 
agents (though constant in time), our preliminary results show 
some patterns analogous to the ones observed here, but with 
much more added complexity. The present investigation lays 
a good foundation on which to base the study of these more 
complicated cases. 

In summary, our results aid in understanding the stability 
of complex coherent structures in swarming systems with 
time delayed communication and in the presence of a noisy 
environment. Although our analytical and numerical results 
were obtained using a model with linear, attractive interactions, 
our analysis gives useful insight for the study of models 
with more general forms of time delayed coupling between 
agents. Our results may prove to be valuable for the control 
of man-made vehicles where actuation and communication 
are delayed, as well as in understanding swarm alignment in 
biological systems. 

Appendix A 
Analysis of the Ring State 

The swarm ring state is obtained when the center of mass 
is stationary. For the solution R =const. to satisfy Eq. (|7]i we 
require 



N 



(29) 



We simplify Eq. (O by taking R =const. and using Eq. 
we obtain 

6r^ = (1 - SPj) Svj - aSvj - ^Srj{t - r). (30) 

We consider the large system size limit — > oo and we drop 
the delayed term. The resulting equations are simply ODEs 
and so the analysis below shows that the ring orbit is not 
dependent on having time delays in the system. Writing Eq. 
( [3Ql l in polar coordinates 5xj ~ pj cosOj and Syj ~ pj cosOj, 
we obtain 



pA = (i-p',-pP'j)pA-2p,9,. 



(31a) 
(31b) 



Equations ( l31al )-( l3Tbl ) have the trivial solution pj = as 
well as a ring solution: 



in which particles move at unit speed, pj9j = ±1. 



(32) 
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Appendix B 
Analysis of the Rotating State 

In the noiseless rotating state, all particles collapse to a 
point, Sr.i = 0, and the equation for the center of mass given 
by Eq. (jTji simplifies considerably to 



(33) 



R = 1 - R - a (R(t) - Kit - t)) , 



We write R = {X, Y) and introduce polar coordinates X 
p cos 9 and Y ~ p sin 6 to obtain 



•^^ p^6'^]p + pO'^ — a[p — pTCOii{w ^Ti], 

(34a) 

p'e = [l- p^ - p^e^^ pf) - 2pe + apr sin(6i - Or), (34b) 

where we've written pr = p{t — r) and 6t = 9{t — t). 
Equations ( |34at -( [34b] i have a circular orbit solution, p ^ pq 
and 9 = ojt + 9o, where 



=a ■ (1 — coswr), 



1 / smwT 



(35a) 
(35b) 



\uj\ V w 

and 6*0 is obtained from the initial conditions. In the main text 
we discuss the behavior of the solutions to Eqs. ( I35ab - (l35bb . 

Appendix C 
Analysis of the Degenerate Rotating State 

When the motion of the whole swarm is initially constrained 
to a line, Eqs. (fTat-(fTb]i dictate that the swarm will remain 
on this line for all times. If the coupling parameter a and/or 
the time delay r are large enough, the resulting motion is a 
degenerate form of the rotating solution in which the swarm 
moves back and forth along a straight line. 

In the case without noise all particles collapse to a point, 
(5r, = 0, and the line along which motion occurs is arbitrary; 
here we use X = Y. The problem reduces to analyzing a 
single delay equation given by 

X = (l- 2X^) X -a {X{t) - X(t - t)) 



(36) 

We find a solution using Fourier analysis. We let 

oo 

X{t)= J2 cne™"*, (37) 



where the coefficients satisfy c„ = c_„* in order to ensure 
that X{t) is a real quantity. Substituting Eq. (|37] | into Eq. 
J, we get for the n-th mode 



2 2 

-n LO On — inUJCn 



+ ^iuf" ^ CiC„iCn-i-m(ni{n - £-m) 

- ac„(l - e-™'^"), (38) 
for n = 0, 1, 2, . . . . The n = equation is 

ceCjnC-e-,nim{i + m) = 0, (39) 



which does not involve cq. Unsurprisingly, cq is undetermined 
since the position of the center of mass may be translated in 
space without modifying the dynamics of the system. 

We now approximate the motion of the center of mass by 
keeping the first three modes. By appropriately choosing the 
time origin, we may take ci to be purely real and positive. In 
contrast, C2 and C3 are complex quantities which we write as 
Ci — |ci|e"^\ for i = 2,3. The equations for the first three 
modes n = 1,2,3 become 



lUJCl 
,3 / Q„3 



+ 2iu}^ (-3c? - 36clcl - 54ci|c3|^ - 24ci|c2| 

— 4aj^C2 = 2iu;c2 

+ 2iw3 (_i08c2|c3p - 36cic;c3 - 24c2|c2p - 



ac2(l 



9c?C3) 
(40a) 



12C2C2) 

(40b) 



— 9aj C3 = 3iajC3 

+ ^iufi (-I8C3C1 - 72c3|c2p - 8IC3IC3P - 12c^ci + c?) 
-ac3(l-e-3'"^). (40c) 

In addition, the condition from Eq.([39]l becomes 

6(C2C5 - C2C3) - Ci(c2 - C2) = 0. (41) 

Separating Eqs. ( I40ab -( l40cl i and Eq. (1411 1 into real and imag- 
inary parts yields a system of seven equations (since the 
real part of Eq. dTIT i is satisfied automatically) for the six 
unknowns: uj, c\, |c2|, 02, |c3| and (/)3. These equations cannot 
be satisfied in general. However, if |c2| = 0, then the equation 
for mode n = 2 [Eq. (I40bb 1 and Eq. (HTb are satisfied 
automatically, leaving four equations: 

-w^ci =icjci + liu? (-3c?54ci|c3p + 9C1C3) 

-aci(l-e"'"^), (42a) 



-9a;^C3 =3iwc3 + liu? (-I8C3C1 - 8IC3IC3I 



(42b) 



- ac3(l - e-^*"^). 

for the four unknowns cj, ci, |c3| and (^j,. Equations ( I42ab - 
(I42bb may be solved numerically and permit one to approxi- 
mate the motion of the center of mass in the form 

X(t) = Y{t) ^ 2ci cosujt + 2|c3| cos(3wt + ^3). (43) 

The frequency of the straight line orbit of the swarm center 
of mass is approximately equal to the frequency of the circular 
orbit in Eq. ( 125 ab . In addition, the amplitude of oscillation of 
the straight line orbit is approximately equal to the radius of 
the circular orbit of Eq. (I25bb divided by a factor of \/6- 
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